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Abstract. - The classical Thomson problem of n charged particles confined to the surface 
of a sphere of radius a is analyzed within the Darwin approximation of electrodynamics. For 
n < n c (a) the ground state corresponds to a hexagonal Wigner crystal with a number of 
topological defects. However, if n > n c (a) the Wigner lattice is unstable with respect to small 
perturbations and the ground state becomes spontaneously magnetized for finite n. 



The Thomson problem, finding the ground state of electrons inside a sphere with a uniform 
neutralizing background, has a time honored position in the history of modern physics [1-7] . 
The original question was posed by J.J. Thomson [8] after his discovery of the electron in 
1897. Thomson conjectured that the knowledge of the positions of the electrons inside the 
atoms is essential to understanding the regularity of the chemical elements in the periodic 
table. At the time, however, proton still had to wait 14 years to be discovered, so in order to 
keep his atom neutral, Thomson was forced to introduce a uniform neutralizing background. 
The model became known as the "plum pudding" atom and the question that needed to be 
answered was: What are the positions of the electrons inside a uniformly (positively) charged 
sphere? Surprisingly, after more than a century this problem still has no general solution. 

If the background charge is made to vanish, the electrostatic energy will be a minimum 
only if all the electrons are located at the surface. This is a general consequence of the Earn- 
shaw theorem [9] which precludes existence of a stable equilibrium with purely electrostatic 
interactions. Curiously, the Coulomb potential is precisely on the border line where this be- 
havior is possible. If instead of 1/r, the electrons would interact by a l/r 1+e potential with 
e > 0, the bulk occupation of the sphere would be energetically favorable for a sufficiently 
large number of electrons [10]. Unfortunately, even the restricted surface Thomson problems 
remains unsolved for an arbitrary number of electrons [11, 12]. 
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In this letter we will show that if the relativistic corrections to the Coulomb law are 
properly taken into account, even our intuitive picture of the ground state as consisting of 
stationary particles located at fixed positions on the surface of a sphere must be abandoned. 
Instead, we find that for sufficiently large electron density, the energy is minimized by the 
particles undergoing a coherent motion and the sphere becomes spontaneously magnetized! 

The starting point for our analysis is the well known Darwin Lagrangian [13-23], which 
takes into account the relativistic corrections to the Coulomb law resulting from the particle 
motion, 

l = — 2 £ - 5 £ 7? + 4? £ 7? ■ «J + • ' • ^) 

Eq. (1) is correct to order v 2 /c 2 . The velocity-dependent correction to the Coulomb energy 
arises from the electromagnetic coupling between the moving particles. Since the Lagrangian 
|JT|| does not contain explicit time dependence, the energy of the system 

i 

is a constant of motion, 

^ = £-™ + 5£7? + ^£7?[«i-i + («*-^-^)]- ( 3 ) 
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The ground state for n electrons on the surface of a sphere of radius a is then determined by 
the minimization of Eq. (J3J . 

We note that if the terms of order 1/c 2 are neglected, we recover the classical formulation 
of the Thomson problem in which the electromagnetic coupling between the electrons is purely 
of the Coulomb form. In this case, the velocity dependent contribution to the Hamiltonian 
is positive or zero, and the ground state corresponds to stationary particles residing at fixed 
positions on the surface of the sphere. For large n, this structure resembles a hexagonal Wigner 
crystal containing some topological defects. In general, however, the 1/c 2 terms can not be 
omitted and a full minimization of Eq. J5J must be performed. To proceed, it is convenient 
to rewrite the energy in adimcnsional form. Defining the reduced displacement and velocity 
as r* — r/a and v* — v/c, the reduced energy becomes 

where r e = q 2 /mc 2 is the classical electron radius and a* = a/r e . 

It is convenient to work in the spherical coordinate system with unit vectors n, 8, 4>. The 
reduced velocity of electron i on the surface of the sphere is then v* — Ovgi + cfrvfo. 

Minimization of E* , Eq. J2J , is performed using a general purpose quasi-Newton method 
where the Hessian update is given by the BFGS formula [24]. Gradients are computed ana- 
lytically. A line search with cubic fit is used with the additional safeguard against evaluations 
beyond light speed. The procedure is highly non-trivial. In fact it is known that already 
for the classical Thomson problem, in the absence of relativistic corrections, there exists an 
exponentially large number of metastable states [25]. Thus, it is quite unlikely that any 
minimization procedure will be able to locate the exact ground state for a large number of 
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electrons. This however, is not of great importance since the metastable states have energies 
very close to that of the exact ground state [10,25]. Performing the minimization of E* we find 
that for reduced surface charge density a* = n/a* 2 such that a* < a* (subcritical region), the 
electrons form a stationary Wigner crystal with some topological defects. Above the critical 
charge density a* > a* (supercritical region), the Wigner crystal, however, becomes unstable 
and a new ground state with moving electrons is formed. In Fig. \\\we show the characteristic 
distribution of particles in this new ground state. The arrows indicate the relative magnitude 
and direction of the particle velocities. The figure shows bands of correlated antiferromagnetic 
velocities that try to adapt to the topology of the sphere. We stress that when the instability 
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Fig. 1 - The velocity field of electrons in the supercritical region, n — 400 and a = 0.7a c . 

occurs, Vi/c -C 1 for all the particles, so that the Darwin Lagrangian remains valid, up to 
quantum corrections. The melting of the Wigner crystal is an example of a classical zero 
temperature phase transition. 

To better understand the nature of the instability of the Wigner lattice, it is convenient 
to rewrite the Darwin energy in a matrix form. Defining a 2n component velocity vector 
V = {«xj«2> ■■■' U n}> Eq. (@J can be rewritten to order 1/c 2 as 

E * ^n+lv T IV + -L V T BV + lY vf + -L Y i- , (5) 

2 2a* 8 ^ 1 2a* ^ r* v ; 

where I is a 2n x 2n identity matrix and D is a position dependent matrix constructed from 
the last term of Eq. Q . The quadratic term in velocity is non-negative if all the eigenvalues 
of the matrix 

A = I+^D, (6) 

are positive. In this case the ground state will have V — and the electrons will be organized 
into a Wigner crystal. On the other hand, as soon as one of the eigenvalues of A becomes 
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negative, the Wigner lattice will lose stability, and a new ground state, with energy below 
that of the Wigner crystal will be established. The phase transition occurs when Xf^ in = 0, 
where A^ in is the minimum eigenvalue of the matrix A. 

It is important to note that the energetic bifurcation of Eq. (JjJJl is simultaneous with the 
dynamical instability of the Wigner lattice. If the Euler-Lagrange equations of motion are 
linearized around the stationary positions of the Wigner lattice, one can show that the Lya- 
punov instability occurs precisely when A loses convexity. Unfortunately, in the supercritical 
region, the equations of motion are differential-algebraic and due to the singularity of A are 
very difficult to integrate numerically [26] . 

To determine the critical charge concentration at which the Wigner crystal loses stability, 
we adopt the following procedure. For a given number of electrons n, the Coulomb energy is 
minimized to determine the positions of all the particles. For purely Coulombic interactions, 
the ground state location of the electrons is independent of the size of the sphere, since a 
scales out of the expression for the electrostatic energy. Once the ground state coordinates 
are known, the eigenvalues A D of the matrix D can be calculated numerically. The criticality 
condition A^ in = is then equivalent to the requirements that X^ in = —a* . In Fig. 0we 
show the result of this procedure. 
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Fig. 2 - The critical radius as a function of n. The solid line is a fit given by Eq. Q 
The points in Fig. H can be very well fitted by 
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Eq. J7J implies existence of a well defined thermodynamic limit for the phase transition, 
limn — > oo, a — > oo and a* — > 5.35. 

To further explore the nature of the ground state for a > a c we define an order parameter 



T; X Vi 



(8) 
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Clearly, fi* is just proportional to the total magnetic moment of the sphere. In the subcritical 
region, the electron velocities are zero and fi* — 0. The value of the magnetic moment in the 
supercritical region is plotted in Fig. |3| We find that if the magnetic moment is scaled with 
n~ 2 / 5 and is plotted as a function of the reduced surface charge concentration a* — er*, all the 
points for different values of a and n fall on the same universal curve, 

g{x) = 0.545 x 1/2 . (9) 

Thus, although locally the orientation of the velocity vectors is antiferromagnetic, globally the 
symmetry is broken and the sphere acquires a net magnetic moment. The magnetic moment 
is sub-extensive and vanishes with a square root singularity as a — > er+. 



0.1 



0.08 




* 

3. 



0.06 



0.04 



0.02 




0.005 



0.01 



0.015 



0.02 



cr. 



Fig. 3 - Magnetic moment in the supercritical region, the data points are the averages over five 
random initial configurations. 



We next proceed to study the extensive property of the electromagnetic energy E* . For 
a < <j c a very accurate expression for the ground state of the Thomson problem [27, 28] can 
be obtained using a simple argument. Consider a uniformly positively charged spherical shell 
on which move n electrons. This problem defines a spherical one-component plasma (SOCP). 
The electrostatic (Coulomb) energy can be written as 



Ft 



SOCP 



Er 



2 2 

q n z 
2a 



q 2 n 2 



(10) 



The first term is the Coulomb energy of interaction between n electrons on the surface of the 
sphere, the second term is the self energy of the positive background charge, and the third term 
is the energy of interaction between n electrons and the background. At zero temperature, 
the classical SOCP will freeze into a hexagonal Wigner crystal (with some topological defects) 
whose energy is 



q 2 n 



Fsqcp = -M^— 



(11) 
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where M is the Madelung constant and d is the characteristic size of the Wigner-Seitz cell, 
ircPn = Ana 2 . Combining Eqs. i|l(J|) an( l we arrive at a very simple expression for the 
Coulomb energy of n electrons on the surface of the sphere 

Eq. (|12|l with M = 1.1046 gives a very accurate fit to the ground state energy of the surface 
Thomson problem with purely Coulomb interactions [10,27]. Note that for a planar OCP [29] 
M = 1.1061, so that the topological defects affect very little the value of the Madelung 
constant. It is also important to notice that although Ec is not extensive, 

is. Therefore, if IS.Ec/n is plotted as a function of n/a 2 for different combinations of n and 
a, all points should fall onto one universal curve, 

f{x) = ~yfc. (14) 

We can now check if this universality also holds for the Thomson problem with the Darwin 
coupling between the particles. That is if 

o 

AE* = E* - 1 - — , (15) 
2a* v ' 

is such that AE* = nf(a*), with f(x) given by Eq. JTU). In Fig\Z\AE*/n is plotted as a 
function of a* for various combinations of n and a* . It is quite surprising that even in the 
supercritical region a* > 5.4, the deviation of E* from the energy of a stationary Wigner 
crystal remains very small. This is in spite of the fact that for a > er c , the velocities of 
individual particles can be quite large. Evidently, the local antiferromagnetic ordering of the 
velocity vectors leads to significant cancellations which diminish the overall contribution of 
the Darwin term to the total energy. 

We have shown that if the relativistic corrections are taken into account, the classical 
Thomson problem of the electrons confined to the surface of a sphere exhibits an electromag- 
netic instability. While for a < <j c , the ground state of electrons is a Wigner crystal with 
some topological defects, for a > cr c , the Wigner lattice is unstable and a small perturbation 
can make the system evolve to a new ground state. This ground state is characterized by a 
local antiferromagnetic order [21-23], but finite net magnetic moment. The surface charge 
concentration at the phase transition has a well defined thermodynamic limit n — > oo, a — ■* oo, 
while a* — > 5.4. This surface charge density, however, is so large that quantum effects must 
be taken into account [30,31]. The relativistic corrections to the Coulomb energy should not, 
therefore, affect the stability of a normal plasma. 
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